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Abstract. We investigate the diffusive properties of energy fluctuations in a one- 
dimensional diatomic chain of hard-point particles interacting through a square- 
well potential. The evolution of initially localized infinitesimal and finite pertur- 
bations is numerically investigated for different density values. All cases belong 
to the same universality class which can be also interpreted as a Levy walk of 
the energy with scaling exponent 7 = 3/5. The zero-pressure limit is nevertheless 
exceptional in that normal diffusion is found in tangent space and yet anomalous 
diffusion with a different rate for perturbations of finite amplitude. The different 
behaviour of the two classes of perturbations is traced back to the "stable chaos" 
type of dynamics exhibited by this model. Finally, the effect of an additional in- 
ternal degree of freedom is investigated, finding that it does not modify the overall 
scenario. 

1 Introduction 

Diffusion in one-dimensional systems of interacting particle may display peculiar features. For 
instance, consider an assembly of impenetrable spheres confined within a narrow channel. If 
the mutual passage of particles is forbidden, the motion of the spheres is necessarily correlated, 
even at long times, because the displacement of a given particle over a long distance necessitates 
a rearrangement of many other particles in the same direction. This is an observed effect, for 
example, in single-filing systems where Pick's law is violated [1]. 

The distinguished signature of those effects is in the long-time behavior of the associated 
correlation functions [2]. As it is known, the latter may display long-time tails leading to ill- 
defined transport coefficients or, more generally, to the breakdown of customary hydrodynamics. 
Indeed, power-law decay of correlations is expected to be a generic feature of one-dimensional 
systems in the presence of conservation laws [3|. One important consequence of such long- 
ranged correlations is that physical properties may significantly depend on the system size and 
that the thermodynamic and infinite-time limits may not commute. For instance, the tagged- 
particle diffusion coefficient in = 1, that is finite for the infinite system is found instead to 
vanish for a finite one (see [4] and references therein) . 

Another important example is the divergence of the thermal conductivity coefficient k, ob- 
served in chains of anharmonic oscillators 0. It is now assessed [B] that in ID k diverges with 
the system's length L 

L" (1) 
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A considerable number of papers have reported anomalies of this type and we refer the reader 
to recent review articles [6l[7] for an account of the literature. Related manifestation of such 
anomalies are the divergence of viscosity in II? lattice gases (§■ and the anomalous scaling of 
Rayleigh and Brillouin peak widths in the hydrodynamic limit detected for the Lennard-Jones 

fluid Enuj. 

For what concerns the actual values of a there are still some controversies. The renormalization- 
group arguments of Ref. [3J give a = 1/3. More recently, a critical revision of the self-consistent 
Mode- Coupling Theory approach pTlfT^ has confirmed this prediction whenever the lowest 
nonlinear term in the theory is cubic in the relevant correlators [13] . On the other hand, when 
the first nonlinearity is quartic, the value of the exponent turns out to be a = 1/2 |141I15|. 
Physically, this latter condition can be achieved when the interaction potential is symmetric 
with respect to the equilibrium position, and corresponds, more generally, to the case of zero 
pressure. 

Many efforts have been made to verify numerically the theoretical predictions. Simulations 
of an anharmonic chain with cubic plus quartic non-linearity yield a = 1/3 [13j and the same 
result has been obtained for a ID gas of elastically colliding particles [TB] with a great accuracy. 
On the other hand, the same analysis carried out for chains with quartic nonlinearities (e.g. 
the FPU-/? model), although exhibits a different scaling behaviour is not able to convincingly 
confirm the above mentioned theoretical predictions T71ll3j. 

Computer simulations of simple toy models is an invaluable way of attacking those problems. 
In particular, one would like to understand the conditions under which anomalies occur and 
to classify possibly universal features. In this sense, the Hard-Point Gas (HPG) [18^, made of 
point particles that interact through elastic collisions, proved very helpful in clarifying some 
aspects of heat conductivity in one-dimensional systems [TM^DUTl 112211^ . A limitation of the 
HPG model is that one cannot tune the pressure and, in particular it is not possible to test 
whether the zero pressure limit indeed corresponds to a different universality class. In order 
to overcome such limitations, in this paper we consider an extension of the HPG model in 
which the distance between neighbouring particles is also limited from above by a square well 
potential. For reference, we will denote this model as the hard-point chain model. Savin and 
collaborators [2l] have already shown that this systems displays anomalous properties. What 
is more important for our purposes is that it an equilibrium state with vanishing pressure and 
we can study how anomalous features change approaching this limit. It should be also recalled 
that in Ref. [24] only the zero-pressure case was considered. 

We indeed confirm the zero-pressure case belongs to a different universality class, although 
it is again difficult to accurately determine the corresponding scaling behaviour. Moreover, by 
studying both infinitesimal and finite perturbations we convincingly find different properties 
which we interpret as a signature of "stable chaos" , i.e. of an irregular behaviour induced by 
discontinuities in phase space, rather than from the exponential amplification of infinitesimal 
perturbations [25]. We must indeed recall that both the HPG and the hard-point chain are not 
chaotic in the usual sense, as the maximum Lyapunov exponent strictly vanishes. 

In section 2, the model and its physical parameters are introduced. In section 3, the diffusion 
properties of both infinitesimal and finite amplitude perturbations are investigated for different 
values of the specific volume (i.e. upon indirectly varying the pressure). The different velocities, 
that can be defined dynamically and thermodynamically, are reviewed in Sec. 4. The properties 
of the energy current are briefly discussed in Sec. 5 to confirm the indirect observation of an 
anomalous behavior at zero pressure. In Sec. 6 the effect of an additional internal degree of 
freedom is analysed, while the last section is devoted to summarizing the results and the open 
problems. 



2 The hard-point chain model 

The model consists of a set of N point-like particles of masses rrii and positions Xi ordered along 
a line. Interactions are restricted only to nearest-neighbour pairs with a square-well potential 
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in the relative distances 

U{xi+i - x^) = 



< Xi^i — Xi < a 
oo otherwise 



The two infinite barriers at Xi-^-i — Xi = and Xi^i ~ Xi ~ a correspond to elastic "collisions", 
the latter ones occurring at a finite distance as if the particles were linked by an inextensible 
and massless string of fixed length a (see Fig.[T]for a pictorial image). The string has no effect on 
the motion unless it reaches its maximal length, when it excerts a restoring force that tends to 
rebounce the particles one against the other. Clearly, the potential ([2]) introduces the physical 
distance a as a parameter of the model. 




Fig. 1. Sketch of a diatomic hard-point chain. Different radii are suggestive of different masses of 
particles that are assumed to be point-like. 



Both types of collisions are described by the same updating rule for the velocities Ui — Xi, 



, rrii - nii+i 

Ui = ■ Ui 



2m 



rrii + rrii+i 



rrii + rrii+i 



I 2mi rrii ~ i^i+i 

Ui+i , Mi+i = : Ui ■ u^+i, (2) 

mi+ nii+i nil + m^+i 



We consider the case of alternating masses, whose ratio is denoted by r, i.e. rrii — ni for 
even i and rrii — , otherwise. This choice is dictated by the need to break the integrability 
of the model that holds for r = 1. The particles are confined inside a "box" of length L. We 
consider both periodic boundary conditions 



Xn+N 



Xn ~\~ L 



and fixed boundary conditions. 



.Ti = Xn — L 



(3) 



(4) 



Accordingly, the specific volume £ = L /N is a state variable to be considered together with the 
specific energy (energy per particle) e — '^^rriiuf /2N . 

As it is well known, the thermodynamics of models like ours can be solved exactly (see e.g. 
|26j). In particular, the partition function and the equation of state in the P — T ensemble is 



Y[ / dpj / dxi / dx2... / dxTv exp[-/3^ — ^ Pxn] 

j^l-J-oo Jo J XI JxN~l i 



(5) 



where pi = rriiUi are the conjugate momenta. Introducing the new variables S^i — a;i+i — Xi we 
find, after some manipulations 



Z = \ 



-N 



Jo 



N 



1 



A/3P 



(1 - exp(-/3Pa)) 



N 



where A is the De Broglie wavelength. The equation of state is given by 



d\nZ 

'Wp) 



N 



1 



PP expiPPa) - 1 



(6) 



(7) 
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This equation is shown in Fig. ([2]) for a = 1. For convenience we report the specific volume £ as 
a function of f3P. Note that for large a values, the equation of state is the same for an HPG (red 
line) i.e. the one of an ideal gas in ID. As a check of our codes (see below) we also computed 
some state parameters by fixing £ (namely the box length) and evaluating the pressure as the 
average exchanged momentum between neighbouring particles. 

1 

UU 
0.8 

0.6 

0.4 

0.2 
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Fig. 2. Equation of state of the hard-point chain: the solid line is the analytic expression of eq[71 points 
are obtained from numerical simulations with e = 1. 

Before discussing the results, let us briefly describe the simulation procedure. Models like 
the one we are dealing with are particularly suitable for numerical computation as they do not 
require integration of nonlinear differential equations. Indeed, the dynamics amounts simply to 
evaluating successive collision times and updating the velocities according to Eqs. The only 
errors are those due to machine round-off. Moreover, the simulation can be made very efficient 
resorting to fast updating algorithms. The fact that the Xi contribute only indirectly to the 
evolution, by determining the collision times together with the conservation of the particle 
ordering along the chain, allows simulating the dynamics with an event driven algorithm that 
exploits the heap structure of the future collision times [21] . 

Without loss of generality, we can fix m = 1 and a — 1. Most runs have been performed for 
a mass ratio r = 2 and for an energy per particle e = 1, which corresponds to = fc^T = 2. 
Moreover, we have chosen to work with two representative values of the specific volume, namely 
£ = 0.2 and £ = 0.5. In the former case, the particles hardly feel the barrier at relative distance 
1, mainly interact only through standard collision. In this regime, the hard-point chain behaves 
like a gas and is indeed practically equivalent to the HPG. On the contrary, for £ — 0.5, the 
two types of collisions are equally probable, so that the pressure is zero. 

3 DifFusion of perturbations 

In the spirit of linear response theory |27| , transport coefficients can be determined by looking 
at the way a statistical system reacts to small perturbations (external forces) . For heat conduc- 
tivity, one can directly follow the evolution of an initially localized energy perturbation Ae(x, t) 
[28j . If the perturbation is weak enough for the system to be locally in thermodynamic equilib- 
rium, a reasonable measure of the local temperature field AT{x,t) is Ae(x,<) = CpAT{x,t), 
where Cp is the heat capacity per unit length at constant pressure. Thus, by measuring the 
spreading of the perturbation field, one can estimate how heat propagates through the system 
[28j . A standard diffusive process is the indication of a normal heat conductivity. 

In general, quantities like Ae(x,<) are determined by all particles contained in a small box 
around the spatial location x. If, like in our ID model, particles cannot cross each other, it 
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is more practical to identify x as the average position it of the i-th particle. To some extent, 
this amounts to interpreting the model as a lattice system. Of course, nothing changes in the 
underlying physics, provided that one consistently adopts the same interpretation throughout. 
Moreover, since we are mainly interested in the long-time and large-distances behaviours this 
should not be very relevant, anyhow. 

With the above remark in mind, we consider a system at equilibrium with an energy eo per 
particle and perturb it by increasing by some preassigned Ae the energy of a subset of adjacent 
particles. Let e(x^t) be the energy profile evolving from such a perturbed initial condition. We 
then ask how the perturbation 

f{x,t) = {e{x,t)-eo) (8) 

behaves in time and space. The angular brackets denote an ensemble average over independent 
trajectories. Because of energy conservation, ^ ■ f{it, t) — Ae at all times, so that / can be 
interpreted as a probability density (provided that it is also positive-defined and normalized). 
At sufficiently long times, one expects that f{x,t) scales as 

f{x,t) = t-'Gix/t'') for \x\ < vt. (9) 

for some finite velocity v. 

A sufficiently general framework where the above scaling behaviour is indeed observed is 
that of Levy walks, where f{x, t) is the probability density of finding a particle whose kinetics 
is defined as follows. The particle moves balistically in between successive "collisions" whose 
time separation is distributed according to a power law, oc t~^~^ , /J. > 0, while its velocity 
is chosen from a symmetric distribution tf'(w). By assuming a (5-like distribution (tf'(w) = {5{v — 
v) + 5[v + v))/2), the propagator P{x, t) (the probability distribution function to find in x at 
time i, a particle initially localized at x = 0) can be written as P{x, t) — Pl{x, t) + t^~^[S{x — 
vt) + 6{x - vt)] where [29 

( t-i/M exp [~{^Jx/t^/^')^] \x\ < t^/" 
PL{x,t) cx \ tx-^'-^ t^^f' < \x\ < it (10) 

[ \x\> vt 

Smoother, choices of the velocity distribution 'F{v) lead to different broader side-peaks, but do 
not affect the shape and the scaling behaviour of the bulk contribution Pl{x, t) which scales as 
predicted in Eq. ([9]) with the exponent 7 = Notice also that the two further parameters 
V and 77 can be scaled out by using proper units of measure for the time and space variables. 
Nevertheless, we shall see in the following that they can carry some useful information when 
approaching a transition point. 

From the evolution of the perturbation profile, it is possible to infer the growth rate a of the 
thermal conductivity (see Eq. ([1])). In fact, in Ref. [SO^, it has been found that a, the scaling 
behaviour of the mean square displacement {<j'^{t) = ^ii^f{x — ii,t) cx t^ and 7 are linked 
by the following relationships. El) 

a = /3-l = 2-i. (11) 
7 

In particular, we see that the case 7=1/2 corresponds to normal diffusion {P — I) and to 
a normal conductivity (a = 0). On the other hand, 7 = 1 corresponds to a ballistic motion 
(/3 = 2) and to a linear divergence of the conductivity {a — 1). 

3.1 Infinitesimal perturbations 

With reference to ID systems, it is convenient to introduce the variable Ui = yfmlui and 
express the perturbed trajectory as = tt^-l-^i, where u^it) is an equilibrium trajectory. Since 

^ Although fi is here different from l/ksT, we use this letter for consistency with previous pub- 
lications. We are anyhow confident that it is always clear from the context which j3 we are referring 
to. 
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in gas-like systems there is only the kinetic contribution to the energy, Eq. ([8]) can be written 
as 

= t) = ((^D^) - = + 2{u%). (12) 

In the limit of infinitesimal perturbations, 8i follows the tangent space dynamics, i.e. the equa- 
tion of motion linearized around a reference trajectory. For the hard-point chain this is for- 
mally equivalent to that in real space, since the dynamical equations are piecewise linear (the 
transformations ([2|) describing the elastic collisions are indeed linear). Accordingly, the en- 
ergy conservation in real space transforms itself into the conservation of the Euclidean norm 
Q = ^if(i,t) = ^i^fii) of ^ generic perturbation. This, in turn, implies that hard-point 
systems are not chaotic in the usual sense of the word (all Lyapunov exponents are equal to 0). 
Moreover, it is reasonable to conjecture (and we have verified numerically) that the dynamics 
in tangent space is uncorrelated with that in real space, i.e. (u^Si) = 0. As a result, the energy 
diffusion can be studied by simply looking at the spreading of infinitesimal perturbations 

f{x^ie,t)^{5^) (13) 

This approach has the advantage of dealing with positive-defined quantities 6f and thus allows 
getting rid of the strong fluctuations affecting positive/negative variables. It is important to 
notice that this is possible precisely because of the lack of standard deterministic chaos: in 
chaotic systems, perturbations would exhibit an exponential growth, which can be cancelled 
only after averaging over an unthinkable number of realizations. 

We have first computed f{x,t) for the hard-point chain with an average interparticle dis- 
tance (specific volume) £ = 0.2, a value sufficiently small to expect a perfect agreement with the 
behaviour of the HPG. In fact, in Fig.O we see that 7 = 0.6 leads to a very good data collapse 
in the central region (the side peaks obviously move away, as their position grows linearly in 
time). 




These results confirm the behaviour of the HPG described in Ref. [16], where it was also 
found a perfect agreement with the scaling function predicted for a Levy walk with fj, — 5/3 
(see Eq. ([To])). In fact, the only appreciable difference with the theoretical prediction, the broad 
secondary peaks exhibited by f{x,t), disappears as soon as the (5-Dirac distribution for V{v) is 
replaced by a Gaussian distribution with a suitable width. 

Upon increasing the specific volume i, the effect of finite-distance collisions (that are the 
novel ingredient of our model) becomes increasingly important. One may thus wonder whether 
the Levy walk scaling ((T0|) still describes the perturbation dynamics. To answer this, in Fig. IH 
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we have plotted the logarithmic derivative 

^^" = dbi^ 

of the maximum height for two values of the specific volume. For i — 0.3 Jeff remains close 
to 0.6, though slightly smaller. As it is implausible that the asymptotic value of 7 depends 
smoothly on any parameter, we attribute the deviations to finite-size corrections. On the other 
hand, the deviations from 0.6 observed for £ = 0.5 are much stronger: the vicinity of jeff to 
1/2 is strongly suggestive of a normal diffusion. 
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Fig. 4. The effective exponent "/eff, Ea. (|14|l . versus time for two different specific volumes: £ = 0.3 
(circles) and 1^0.5 (diamonds), N = 2045 



In order to further clarify the dependence on the specific volume, we have plotted in Fig. [5] 
the profiles corresponding to three different specific volumes. The x axis is scaled so as to yield 
a unit variance. The good overlap among the various profiles confirms that all cases where the 
pressure differs from zero exhibit the same universal behavior (the different position of the side 
peaks reflects simply the different timing and the different ballistic velocities in the three cases). 




-5 5 ^ 



Fig. 5. The rescaled profiles for £ = 0.2 (at time t = 160), £ = 0.3 {t = 280), and £ = 0.4 {t = 580). In 
each case x and / are scaled to yield a unit variance and unit area. 
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This is to be complemented by the observation that the perturbation profile exhibits, in 
the zero-pressure limit a clearly different behaviour (see Fig. [6] where the profiles are scaled 
according to 7 = 1/2). No side peaks are present in this case and the bulk appears to convergence 
towards a Gaussian profile (see the dashed curve, added for reference). 

Altogether, all these results suggest that the zero-pressure case is singular and belongs to a 
different universality class as expected from the mode-coupling theory. This is further confirmed 
by the behaviour of the parameter rj (see Eq. ()10p ) . A convenient way to define rj is by scaling x 
in such a way that the variance of f{x/f )f is equal to 1. As a result, we find that 'q{t) diverges, 
while approaching £ = 1/2. A fit of the few available points is consistent with an inverse square 
root divergence, but many factors prevent considering this prediction too seriously. 



10 




Fig. 6. Rescaled perturbation profiles for I = 0.5, N = 2045, at t = 50, 100, 200, 400. and 750. 
The profiles have been obtained by averaging over 10* realizations. The dashed line corresponds to a 
Gaussian. 



Although this result is a clear confirmation that the zero-pressure case corresponds to a 
different universality class, it poses a problem, since normal diffusion should correspond to 
normal conductivity. This prediction is totally in contrast with both the numerical results for 
the FPU-/3 system and the solution of the mode-coupling equations which instead suggest a 
stronger divergence for the conductivity. 

3.2 Finite perturbations 

In order to understand the puzzling behavior of infinitesimal perturbations, it is worth exploring 
the dynamics of finite rather than infinitesimal perturbations. In fact, the signature of "stable 
chaos" , a type of irregular behaviour found in linearly stable coupled map lattices, is precisely 
the different propagation rate for finite and infinitesimal perturbations 31j. Moreover, since it 
was already argued in Ref. [H] that the HPG might be a Hamiltonian model exhibiting stable 
chaos, it is not unreasonable to expect a different behaviour for the two classes of perturbations. 

With reference to Eq. (fT2|) . we follow pairs of trajectories {u^ {t) , {t)) , where the first 
element is a typical unperturbed equilibrium configuration and the second element is obtained 
from the first one by adding a given amount of energy in the three sites around the origin. 
We then study the evolution of the perturbation by computing the ensemble average of the 
instantaneous local distance A{x,t) = ((uf (i) — u^{t)^) and the ensemble average of the energy 
difference Ae{x,t) = {{uf{t) f} - ((u°(i))'>- 

From the results reported in Fig. [71 the central peaks show a good overlap for 7 — 0.55. 
Considering the shorter time scales with respect to those ones achieved in the corresponding 
tangent-space simulations, the deviation from 0.6 cannot be taken as an indication of a different 
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scaling behaviour. A substantial difference is instead found by comparing the velocity of the 
side peaks which is definitely larger in this second context (see the next section for a detailed 
discussion). The too short time scales prevent a quantitative comparison with the predictions 
of the Levy walk model: in fact the central peak should decay by another order of magnitude 
to appreciate deviations from the Gaussian behaviour. 




Fig. 7. Energy perturbation profiles for an energy density e = 1, at t — 5, 10, 20, and 35 (A'^ — 2045), 
for £ — 0.2. The profiles have been obtained by averaging over 10* realizations. 



In Fig. [8] we plot the amplitude of the perturbation itself. The propagation velocity is 
estimated by determining the half-height position of the front (see the horizontal dashed line) 
at different times. The velocity c/a turns out to coincide with that of the side peaks in the 
set-up of Fig. [7l 
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X 

Fig. 8. Propagation of finite amplitude perturbations for TV = 2045, energy density e — 1 and £ — 0.2. 
The profiles have been obtained by averaging over lO'' realizations. Different curves are at times t = 5, 
10, 20 and 35 respectively. 



Finally, we have studied the case £ = 0.5. As it can be seen in Fig. [9l in the available 
time range we do not appreciate any difference with respect to the previous non-zero pressure 
case. Moreover, in agreement with tangent space results, the side peaks are basically absent. 




10 



Will be inserted by the editor 



Altogether, fluctuations in real space are so strong that it is almost impossible to extract any 
useful information, apart from the velocity of the perturbations themselves, that is the only 
observable unambiguously different with respect to that of infinitesimal perturbations. 




4 Ballistic phenomena 



In the previous section we have encountered two different velocities for the propagation of 
infinitesimal and finite perturbations. From a thermodynamic point of view, there exist two 
further velocities: the adiabatic Cad and the isothermal Cis sound velocity, which are defined 
respectively as 

'"^ \dpJs KdpJr 
where p is the mass density. They are related by 



(15) 



'-'ad 



Cp 
Cv 



(16) 



where Cp and Cv are the constant-pressure and constant-volume specific heats. 

From the equation of state ([5]), we can determine both velocities for the hard-point chain. 
They are reported in Fig. [10] versus the specific volume. For ^ ^ 0, the hard-point chain reduces 
to the HPG model. In this limit, a simple expression can be obtained for the two velocities. 
The equation of state reduces to 



PL = UbNT = 



2p kgT 

1 + r m 



where p - 



m{l + r)L/2N. From Eq. (fT5|) we obtain, 

2 2kBT 

(l+r)m 



(17) 



(18) 



because CpjCy = 3 in the HPG. By inserting the parameters used in the the numerical 
simulations we obtain = 1.154 and Cad = 2 for r = 2, and = 1, Cad — 1.722 for r = 3. 

We can now compare all such velocities. From the position of the side peaks for different 
specific volumes, we obtain the tangent-space velocity Cta (circles in Fig.lTU]). In the HPG limit, 
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Cta coincides with Cis (we have checked that this holds true also for different mass ratios) . How- 
ever, upon increasing the specific volume Cis increases, while Cta decreases and vanishes when i 
approaches 0.5. Besides confirming the peculiarity of the zero-pressure case, these simulations 
reveal that cta differs from both thermodynamical velocities. 

Furthermore, we have determined the velocity c/a of finite-amplitude perturbations (dia- 
monds in Fig. [TO|l . finding that it is close to, though slightly larger than Cad ■ The difference 
is larger than the statistical fluctuations, but we cannot exclude that this is due to systematic 
deviations that vanish at longer times. We have verified that Cfa agrees also with the velocity 
estimated from the dispersion of long-wavelength Fourier modesH 



4 1 1 1 1 1 1 1 1 1 1 r 

C 



3- 




Fig. 10. Thermodynamic velocities versus the specific volume I = L/N for the hard point chain. Solid 
and dashes curves correspond to the adiabatic and isothermal velocity, respectively. The solid circles 
correspond to the propagation velocity as determined from the behaviour of infinitesimal perturbations, 
while the diamonds correspond to the velocity of finite amplitude perturbations. 



From the dynamical point of view, the most important observation is however that c/o is 
definitely larger than Cta. This inequality is the typical signature of "stable chaos". In order 
to strenghten this conclusion we show in Fig. [TT] another typical feature of stable chaos, i.e. 
that self-sustained irregular behaviour is accompanied by the presence of discontinuities in 
phase-space dynamics. By discontinuiy, we mean a codimension one manifold across which an 
infinitesimal difference in the initial conditions is amplified to a finite value in a finite time 
0. In Fig. [TTl we show that this indeed occurs in the HPG in the vicinity of a three-body 
collision. By shifting the central particle to the right, one discontinuously passes from the 
sequence of collisions depicted in panel a to the sequence depicted in panel b. Since the effects 
of the collisions do not commute, one has a discontinuous change in the outcoming velocities, 
while passing through the configuration which yields the three body collision. In the vicinity of 
such discontinuities, small but finite perturbations are suddenly amplified and this provides a 
mechanism for sustaining irregular behavior in the absence of an exponential amplification of 
infinitesimal differences. 



^ The agreement claimed in Ref. [16] between the tangent space and adiabatic velocity was the result 
of an error in the normalization. 

^ Actually, as shown in Ref. [32], the condition is not so strict and strong localized nonlinearities may 
suffice. 
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Fig. 11. The evolution of three neighbouring particles in the vicinity of a three body collisions. 



5 The energy current 

In view of the difRculties encountered in extracting the scahng behavior of thermal conductivity 
from direct energy diffusion simulation, we now turn our attention to the fluctuations of the 
total heat flux. Our past experience with other ID models in fact suggests that this observable 
is much more reliable since the fluctuations due to the spatial dependence are integrated out. 

In order to properly treat hard-point chains, it is necessary to implement the general mi- 
croscopic definition of the heat flux [6] 

ji = ^{^i+i - Xi){ui+i + Ui)F{xi+i - Xi) + im^if = jf^ + jf ^ (19) 

where F denotes the force exerted by the i-th particle on its right neighbor. In the HPG, only 
the second term matter, since interactions occur only when the relative distance between the 
neighbouring particles Xi+i — vanishes. On the contrary, in the hard-point chain, "collisions" 
occur also at the distance a. However, since the force is singular, one cannot implement the 
definition ^9]) as it stands. By defining the force between two particles as the momentum 
difference induced by a collision, jp-* can be written as the kinetic energy variation times the 

actual distance a, i.e., — ami{u^ — u^)/2 divided by a suitable time-interval At. In order 
to get rid of the microscopic fluctuations, it is necessary to consider a sufficiently long At, 
so as to include a large number of collisions. Since the number of collisions is proportional to 
the system size, it is only in long systems that fluctuations can be removed without spoiling 
the slow dynamics of the heat flux. By flxing At = 1 in our arbitrary units, it is sufficient to 
simulate a few hundreds of particles. 

Besides computing the power spectrum S{v) of the total ffux, we have determined the power 
spectrum of the two separate contributions j*^^^ and j^^^^ for i = 0.5. The j'-^-' contribution, which 
exhibits v~'^/'^ power-law divergence in the HPG, here appears to saturate (see the lowest solid 
curve in Fig. [T^ . From the comparison with the spectrum obtained for a shorter chain length 
(see the dashed curve), it is clear that the saturation is not a finite-size effect. Accordingly, 
one might be tempted to conclude that the zero-pressure case is characterized by a normal 
conductivity as also suggested by the simulation performed in tangent space. However, the 
spectrum of the other contribution shows that this is not such conclusion is incorrect. 
In fact, the middle solid curve exhibits a low- frequency divergence which is terminated by a 
finite-size saturation (as revelead by looking at the spectra corresponding to two different chain 
lengths). A power-law fit within the central frequency range (i.e. away from the saturation 
region) yields an exponent around 0.45 (see the straigh segment, shifted for the sake of clarity). 
This value is not only deflnetely larger than the 1/3 value predicted for nonzero pressure, but 
it is also fairly close to the results found for the FPU-/3 model [17]. Although we cannot claim 
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Fig. 12. Power spectrum of the heat flux for energy density e — 1 and specific volume £ = 0.5. Solid 
and dashed lines refer to A'' = 2045 and A'' — 1021 respectively. From bottom to top: the contribution 
due to j'^-*, j'^-* (see Eq. (|19p . and their sum. Power spectrum have been obtained by averaging over 
10* realizations. 



a quantitative confirmation of the 1/2 prediction of mode-coupling theory, the conjecture that 
zero-pressure systems do belong to a different universality class is further strengthened. 

On the other hand, one should not forget that the heat flux definition requires summing the 
two contributions we have separately analysed. By doing so, not unexpectedly, the resulting 
power spectrum exhibits an intermediate behaviour (see the top curve in Fig. I12p . It is reason- 
able to conjecture that the resulting slower convergence is a finite-size correction due to the 
strong constant background due to the j*^^^ contribution, but it is also honest to admit that 
these results altogether confirm the difficulty of extracting reliable estimates in zero-pressure 
systems. 

6 Adding internal degrees of freedom 

For the sake of completeness, we have explored another variant of the HPG that is obtained 
by introducing an additional degree of freedom to each particle in the original model. As it 
is sketched in Fig. [T21 particles are replaced by infinitesimally thin disks, which can rotate 
along the x axis and thereby possess an angular velocity as well as an angular momentum 
that we denote with Q. Like in the original HPG model, collisions do occur whenever the 
longitudinal position of two neighbouring particles coincides and the dynamics is determined 
by the corresponding collision rules. Besides linear momentum Pi and energy Ei, there is another 
conserved quantity, namely angular momentum Zi . For the sake of simplicity, we consider an 
equal-disk gas, while the units of measure are chosen in such a way that both the mass and the 
momentum of inertia are set equal to 1. Altogether, the quantities 



are left invariant by the collision. One can easily verify that the following transformations satisfy 
the above constraints. 



Pi ^ Pi+ Pt+i 

Zi = Ci + Ci+1 



(20) 
(21) 



Pt+i 



c 



Pi sin^ (f) + pi+i cos^ 4) + iCi+i 
Pi cos^ (j) + pi+i sin^ (f) - (Ci+i 
Ct cos^ (f) + Ci+1 sin^ (/) + (pi+i 
Ci sin^ (f> + Ci+i cos^ (j) - {pi+i 



Ci ) sin 4> cos 4> 
(i) sin cos 
Pi) sin cos 
Pi) sin cos 



(22) 
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where (f) is a free parameter. One might determine a deterministic updating rule for by 
assuming a specific interaction Hamiltonian during the collision. For the sake of simplicity, we 
have preferred to randomly extract (jj according to a flat distribution restricted to the angles 
which guarantee that the outgoing velocities change direction (to avoid particle crossing). 

Although this model has been developed by assuming that Ci is an angular momentum, it 
equally applies to the case where Q is a linear momentum, like for a chain embedded in a planar 
geometry with transverse periodic boundary conditions. Actually, this latter interpretation 
would make the model almost identical to that one studied in [33j . 




Fig. 13. Sketch of the hard disk model 



Like before, we have studied the evolution of an energy packet in tangent space, /(x,t) = 
5p1{t) + 5C,f{t), having the care of referring to the same sequence of </> values in the real as well 
as in tangent space. 

The quantitative result of numerical studies are very similar to the diatomic HPG [16j and 
to the hard-point chain model with I = 0.2. More precisely, the evolution in the tangent space 
is in a good agreement with the Levy-walk model with the same scaling exponent 7 — 0.6 
{Yig.Mi- 

Numerical value of the velocity propagation is consistent with the isothermal sound velocity 
of the hard-point chain model with t = 0.2 and equal masses. 

7 Conclusions and open problems 

In this paper we have seen that the hard-point chain model allows extending the same ap- 
proaches developed for the hard-point gas to a system, where the pressure can be tuned from 
positive to negative values. We have restricted our analysis to nonnegative pressure values, 
because the evolution for negative pressures is symmetrically related to that for positive values 
(see, e.g., the equation of state as plotted in Fig. [2]). Although, the various numerical techniques 
have confirmed the difficulties in extracting reliable estimates of scaling rate, we have been able 
to convincingly show that the zero-pressure case belongs to a different universality class. The 
most compelling evidence is arising from the tangent space evolution: (i) the scaling exponent 
7 is close to 1/2 instead of 3/5 as for the hard point gas; (ii) the velocity Cta goes to zero. 

Another general observation follows from the different propagation velocity exhibited by 
finite and infinitesimal perturbations. We interpret this fact as an indication that the behavior 
of hard-point chains is an instance of "stable chaos" , i.e. an irregular dynamical regime sustained 
by discontinuities in phase-space rather than by the exponential amplification. 

Besides these two general comments, all conclusions about the scaling behaviour require 
further checks and a more careful analysis (if doable at all). For instance, although it is quite 
reasonable to assume that the scaling rate 7 is always equal to 3/5 whenever the pressure is 
different from 0, it should be verified that the systematic deviations observed for the specific 
volume I = 0.3 decay when going to longer times. 

On the other hand, the difficulties encountered while performing simulations with finite 
amplitude perturbations indicate that there is little room to make progress, by adopting this 
approach. In this sense, it is perhaps a better idea looking directly at the behaviour of the 
correlation function of the energy density [34| , 



Ce(x, t) = {Se{y, r)(5e(x + y,t + r)) 



(23) 
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Fig. 14. Reseated perturbation profiles, for the evolution in the tangent space, at t — 80, 160, 240 for 
the generalized HPG, Eas. (|22|l . The profiles have been obtained by averaging over 10* realizations. In 
the inset, the profile at t = 160 (solid line) is compared with the propagators of a Levy walk for an 
exponent jj, = 5/3 with a Gaussian distribution with average equal to 1 and rins 0.0347) (dashed line). 



where the angular brackets denote a spatial as well temporal average. At t = 0, the corre- 
lation function Ce{x,0) is a (5 function in space. Moreover, in the microcanonical ensemble, 
energy conservation implies that the area / dxCe{x,t) is constant at all times. By therefore 
assuming that Ce{x,t) is normalizet to have unit area, its behaviour is formally equivalent to 
that a diffusing probability distribution. This observation has led Zhao to determine the scal- 
ing behaviour of the heat conductivity from the growth rate of the variance of Ce{x,t) [34] . 
As the determination of the variance is troubled by the fluctuating tails, it is preferrable to 
proceed like in the previous sections, by looking at the decay of the maximum Ce(0, t) that is 
statistically more reliable. For the sake of providing a concrete example of such an approach, 
we performed a microcanonical simulations of an FPU-/3 chain in a lattice of 2048 particles. 
The correlations Ce{x,t) for different times, are reported in Fig. 1151 They have been reseated 
by choosing 7 — 0.642, a value that has been determined from a best fit of the maximum of 
the correlation function. From the Levy- walk machinery, it follows that a = 2 — 1/^ k, 0.44. 
This value is significantly different from 1/3 and definitely confirms that zero- pressure systems 
do belong to a different universality class. On the other hand 0.44 is still not so close to the 
mode-coupling prediction (1/2) and thus confirms once more the presence of strong finite-size 
corrections which certainly need be understood before making conclusive claims. 
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